Packages

library(dplyr)
library(dygraphs)
library(forecast)
library(h2o)
library(imputeTS)
library(lubridate)
library(plotly)
library(readxl)
library(TSstudio)

Univariate Time Series Data for U.S. Beef Prices

Data set: U.S. average price of 100% ground beef in dollars per pound Period: January 1, 1984, to March 1, 2022 Seasonally adjusted: No Source: Economic Research data from the Federal Reserve Bank of St. Louis (FRED) Series: APU0000703112 Web address: https://fred.stlouisfed.org/series/APU0000703112

This univariate time series data set contains the average monthly values for the price of 100% ground beef in the United States, measured in dollars per pound, from January 1st, 1984, to March 1st, 2022, for a total of 459 months. The missing value for October 1st, 2012, was imputed using linear interpolation.

Exploratory Data Analysis

Plot of observed values

Decomposition of time series data

## [1] "decomposed.ts"
##          Length Class  Mode     
## x        459    ts     numeric  
## seasonal 459    ts     numeric  
## trend    459    ts     numeric  
## random   459    ts     numeric  
## figure    12    -none- numeric  
## type       1    -none- character

The time series has a growing trend with an embedded cycle, which are both apparent in the observed series. The most recent cycle started just before 2010, near the end of the Great Recession that began in 2008. The seasonal component is not apparent in the observed series. The impact of the COVID-19 pandemic from 2020 to 2022 is conspicuous in both the observed series and the random component. The time series plot can be decomposed to show the trend (including cycle), seasonal, and random components.

Seasonality analysis

The heatmap shows evidence of cyclic behavior (across the vertical bars), but not seasonal behavior (across the horizontal bars). All four seasonality plots lack evidence for seasonal behavior in the time series based on the following behavior:

  • Horizontal lines in the standard plot
  • Rope appearance in the cycle plot
  • Horizontal pattern across the box plots
  • Circular spiral pattern in the polar plot

Correlation analysis

## 
##  Box-Ljung test
## 
## data:  bf_ts
## X-squared = 453.32, df = 1, p-value < 2.2e-16

The correlation of the series with its lags is decaying gradually over time, with no apparent seasonal component.

The lack of seasonality makes sense given that beef is a food eaten year-round in the United States.

Forecasting Strategies

Linear Regression

The time series data set was partitioned into a training set consisting of the values of the first 447 months, and a test set consisting of the last 12 months.

## 
## Call:
## tslm(formula = bf_train ~ season + trend + I(trend^2))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.5603 -0.1501 -0.0601  0.1336  0.8855 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.402e+00  5.216e-02  26.874  < 2e-16 ***
## season2      3.606e-03  5.711e-02   0.063    0.950    
## season3     -2.804e-03  5.711e-02  -0.049    0.961    
## season4      8.675e-03  5.750e-02   0.151    0.880    
## season5      5.713e-03  5.750e-02   0.099    0.921    
## season6      1.722e-02  5.750e-02   0.300    0.765    
## season7     -3.282e-03  5.750e-02  -0.057    0.955    
## season8      4.739e-03  5.750e-02   0.082    0.934    
## season9     -2.877e-03  5.750e-02  -0.050    0.960    
## season10    -1.752e-02  5.750e-02  -0.305    0.761    
## season11    -1.962e-03  5.750e-02  -0.034    0.973    
## season12    -1.600e-02  5.750e-02  -0.278    0.781    
## trend       -2.544e-03  3.659e-04  -6.952 1.33e-11 ***
## I(trend^2)   2.082e-05  7.909e-07  26.325  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2489 on 433 degrees of freedom
## Multiple R-squared:  0.9349, Adjusted R-squared:  0.933 
## F-statistic: 478.4 on 13 and 433 DF,  p-value: < 2.2e-16
##                         ME      RMSE       MAE        MPE     MAPE      MASE
## Training set -4.374877e-17 0.2450012 0.1840790 -0.9673158 8.353438 1.2475944
## Test set     -3.864324e-02 0.1766516 0.1341164 -1.0490354 3.074494 0.9089732
##                   ACF1 Theil's U
## Training set 0.9676458        NA
## Test set     0.6193823  1.281222

## 
##  Ljung-Box test
## 
## data:  Residuals from Linear regression model
## Q* = 4848.4, df = 10, p-value < 2.2e-16
## 
## Model df: 14.   Total lags used: 24

The model coefficients for the intercept and the trend components are statistically significant at a level of less than 0.001. None of the coefficients for the seasonal components are statistically significant.

The MAPE is 8.35% on the training set and 3.07% on the test set. This imbalance suggests that the predictors may be overfitting the model, consistent with the high adjusted R-squared value of 0.93. Alternatively, the imbalance may be due to the dates of the test set occurring during the global pandemic, which is an event that is unprecedented in the time period of the data set from 1984 to the present, resulting in outliers.

The Ljung-Box test has a statistically-significant p-value, rejecting the null hypothesis that the random component is white noise. This is confirmed by visual analysis of the residuals, which show significant correlation in the model between the series and its lags. This means that the model does not capture a majority of the variation patterns of the series. Therefore, it is not a valid model for consideration. However, we will use its MAPE score of 3.07% as a benchmark to evaluate the performance of the other models that we will train.

Exponential Smoothing Models

Holt-Winters model

The time series data set was partitioned into a training set consisting of the values of the first 447 months, and a test set consisting of the last 12 months.

## Holt-Winters exponential smoothing with trend and additive seasonal component.
## 
## Call:
## HoltWinters(x = bf_train)
## 
## Smoothing parameters:
##  alpha: 0.9402085
##  beta : 0.008442693
##  gamma: 1
## 
## Coefficients:
##             [,1]
## a    4.005208800
## b    0.007774959
## s1   0.040817459
## s2   0.037682033
## s3   0.036451698
## s4  -0.033699122
## s5  -0.017816973
## s6  -0.019629350
## s7  -0.029568808
## s8  -0.010102702
## s9  -0.011698607
## s10  0.028851510
## s11  0.036848593
## s12  0.036791200
##                       ME       RMSE        MAE       MPE     MAPE     MASE
## Training set 0.003649732 0.06282926 0.04206687 0.1870582 1.878267 0.285108
## Test set     0.427510057 0.47348183 0.42751006 9.3209263 9.320926 2.897448
##                    ACF1 Theil's U
## Training set 0.04762175        NA
## Test set     0.64714940  3.944431

The Holt-Winters model is mainly learning from the level and seasonal update (with \(\alpha=0.94\) and \(\gamma=1\)). On the other hand, there is essentially no learning from the trend value (\(\beta=0.008\)). The accuracy metrics of the model are imbalanced, with an MAPE of 1.9% in the training set and 9.32% in the testing set. This makes sense given the global pandemic and overlapping period of acute inflation occurring from early 2020 to the present, which have no precedent since the beginning of the training set in 1984. This can be seen in the plot of model performance. While the plot shows a good fit to the peak in 2020, most of the error comes from underestimating the following peak in 2021.

Residual analysis suggests that the residuals are white noise, so we can conclude that this is a valid forecasting model.

We will next train the Holt-Winters model using a grid search. We will start with a shallow search with larger increments for the tuning parameters. This will narrow the search area for a deeper search of the tuning parameters. The shallow search will consist of backtesting on the training data set, with an expanding window of six different periods spaced six months apart. Each of the three tuning parameters will be initialized to a range of 0 to 1, with an increment of 0.1.

## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo 
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo 
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
##    alpha beta gamma         1         2         3         4        5        6
## 1    0.2  0.2   0.3 1.0556192 0.7992706 0.9011382 0.8689359 5.289649 4.418293
## 2    0.2  0.2   0.2 1.0074986 0.7718232 1.0665280 0.8438324 5.163495 4.515609
## 3    0.2  0.2   0.4 1.2023497 0.8443070 0.8916600 0.8712494 5.346879 4.403994
## 4    0.2  0.2   0.1 0.9426384 0.7665015 1.3735227 0.9103540 4.932057 4.802884
## 5    0.2  0.2   0.5 1.4010731 0.9159756 0.9363370 0.8641104 5.379831 4.423941
## 6    0.2  0.2   0.6 1.6855986 0.9731892 1.0153964 0.8563282 5.434246 4.443318
## 7    0.3  0.2   0.4 2.2531333 1.5547452 0.8915947 0.9651059 5.561521 4.081708
## 8    0.2  0.2   0.7 2.1472002 1.2228704 1.0708194 0.8680537 5.556067 4.445365
## 9    0.3  0.2   0.5 2.1400805 1.5487270 0.9467306 0.9295308 5.667935 4.092696
## 10   0.3  0.2   0.6 2.0456472 1.6412760 1.0050017 0.9133077 5.667517 4.168704
##        mean
## 1  2.222151
## 2  2.228131
## 3  2.260073
## 4  2.287993
## 5  2.320211
## 6  2.401346
## 7  2.551301
## 8  2.551729
## 9  2.554283
## 10 2.573575

The table displays the results of the shallow grid search. The models are sorted from best to worst, according to the combination of tuning parameters having the lowest mean error rates. The optimal range of \(\alpha\) varies between 0.2 and 0.3, \(\beta\) is constant at 0.2, and \(\gamma\) is between 0.1 and 0.7. This will help us narrow the ranges of parameter values for a deeper grid search.

## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo 
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo 
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
##    alpha beta gamma         1         2         3         4        5        6
## 1   0.22 0.18  0.24 0.8681091 0.7160404 0.9245224 0.8385599 5.155048 4.441463
## 2   0.22 0.18  0.23 0.8660223 0.7144971 0.9402757 0.8345944 5.136618 4.453896
## 3   0.22 0.18  0.25 0.8735209 0.7202531 0.9091592 0.8421694 5.172506 4.430117
## 4   0.22 0.18  0.22 0.8633114 0.7131452 0.9564106 0.8302542 5.117214 4.467424
## 5   0.22 0.18  0.21 0.8599205 0.7119471 0.9729181 0.8255214 5.096840 4.482350
## 6   0.21 0.18  0.27 0.8537151 0.7032510 0.9320996 0.8312069 5.198051 4.436071
## 7   0.22 0.18  0.26 0.8796032 0.7275358 0.8963372 0.8454425 5.189002 4.419839
## 8   0.21 0.18  0.26 0.8641353 0.6892276 0.9474546 0.8274143 5.184209 4.445618
## 9   0.22 0.18  0.20 0.8557948 0.7108650 0.9897899 0.8218965 5.075504 4.505572
## 10  0.21 0.18  0.28 0.8521247 0.7167731 0.9171133 0.8346831 5.211033 4.427773
##        mean
## 1  2.157290
## 2  2.157651
## 3  2.157954
## 4  2.157960
## 5  2.158249
## 6  2.159066
## 7  2.159627
## 8  2.159676
## 9  2.159904
## 10 2.159917

The error range of the top 10 models has dropped slightly compared to the shallow search. The next step is to retrain the HW model using the optimal values of the smoothing parameters from the deep grid search.

## Holt-Winters exponential smoothing with trend and additive seasonal component.
## 
## Call:
## HoltWinters(x = bf_train, alpha = deep_grid$alpha, beta = deep_grid$beta,     gamma = deep_grid$gamma)
## 
## Smoothing parameters:
##  alpha: 0.22
##  beta : 0.18
##  gamma: 0.24
## 
## Coefficients:
##             [,1]
## a    4.086383981
## b   -0.010245597
## s1   0.034766056
## s2   0.090230495
## s3   0.139312464
## s4   0.016836530
## s5  -0.006100014
## s6  -0.033286042
## s7  -0.059159900
## s8  -0.052489039
## s9  -0.069468318
## s10 -0.040308926
## s11 -0.033346079
## s12 -0.012600020
##                         ME       RMSE        MAE        MPE      MAPE      MASE
## Training set -0.0003145299 0.09035163 0.06018836  0.0350331  2.516007 0.4079263
## Test set      0.4735134660 0.55655997 0.48518475 10.2513374 10.536007 3.2883375
##                   ACF1 Theil's U
## Training set 0.6989335        NA
## Test set     0.7380257  4.607591

As you can see from the plot of fitted and forecasted values, the HW model obtained from the grid search underestimated the 2021 peak, with an MAPE score of 10.5% on the test set. Correlation analysis suggests that the random component is not white noise, meaning that this is not a valid model for consideration.

Forecasting with ARIMA Models

Transforming a non-stationary series into a stationary series

The log transformation with first-order differencing did the best job of transforming the series to a stationary state and stabilizing the series variation. We will next use this sequence of transformations to manually fit an ARIMA model.

Fitting an ARIMA model with a manual process

Autoregressive Integrated Moving Average (ARIMA)

## 
## Call:
## arima(x = log(bf_ts), order = c(1, 1, 0))
## 
## Coefficients:
##           ar1
##       -0.1019
## s.e.   0.0466
## 
## sigma^2 estimated as 0.0005695:  log likelihood = 1060.93,  aic = -2117.86
## 
## Training set error measures:
##                       ME      RMSE        MAE       MPE     MAPE      MASE
## Training set 0.003126903 0.0238378 0.01758745 0.3047038 3.311766 0.9898104
##                     ACF1
## Training set -0.01944971

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(1,1,0)
## Q* = 19.599, df = 23, p-value = 0.666
## 
## Model df: 1.   Total lags used: 24
## 
## Call:
## arima(x = log(bf_ts), order = c(0, 1, 1))
## 
## Coefficients:
##           ma1
##       -0.1043
## s.e.   0.0464
## 
## sigma^2 estimated as 0.0005693:  log likelihood = 1061,  aic = -2118
## 
## Training set error measures:
##                       ME       RMSE        MAE       MPE     MAPE     MASE
## Training set 0.003166816 0.02383408 0.01758675 0.3090378 3.310219 0.989771
##                     ACF1
## Training set -0.01698237

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(0,1,1)
## Q* = 19.008, df = 23, p-value = 0.7008
## 
## Model df: 1.   Total lags used: 24
## 
## Call:
## arima(x = log(bf_ts), order = c(1, 1, 1))
## 
## Coefficients:
##          ar1      ma1
##       0.0181  -0.1220
## s.e.  0.3041   0.2996
## 
## sigma^2 estimated as 0.0005693:  log likelihood = 1061,  aic = -2116.01
## 
## Training set error measures:
##                       ME       RMSE        MAE       MPE     MAPE      MASE
## Training set 0.003172063 0.02383404 0.01758691 0.3096127 3.310142 0.9897797
##                     ACF1
## Training set -0.01731726

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(1,1,1)
## Q* = 18.951, df = 22, p-value = 0.6483
## 
## Model df: 2.   Total lags used: 24

A log transformation was applied to the series, followed by first-order differencing. The transformed series cuts off after the first lag in the plots of both the autocorrelation function (ACF) and the partial autocorrelation function (PACF). There is no apparent seasonality pattern. The lags do not appear to tail off in either plot. Therefore, the corresponding models that were trained are ARIMA(1,1,0), ARIMA(0,1,1), and ARIMA(1,1,1).

All three models have nearly identical performance in terms of error metrics and behavior of the residuals. The ARIMA(1,1,1) had a MAPE score of 3.31%, being the lowest of the three models by one thousandth of one percent. Furthermore, the ARIMA(1,1,1) had the most statistically significant coefficients of the three models. The plot of residuals from the ARIMA(1,1,1) model show a pattern of random oscillation around zero with stable variation. There are a few outliers, with the most obvious one corresponding to the global pandemic from 2020 to 2022. Otherwise, there is apparent nonrandom pattern in the residuals. The ACF plot confirms the lack of significant autocorrelation with lagged values. The density plot shows that the errors are normally distributed. The Ljung-Box test failed to reject the null hypothesis that the lags are not correlated. Therefore, we conclude that the random component of the model is white noise, meaning that ARIMA(1,1,1) is a valid forecasting model for this data series.

Fitting an ARIMA model with an automated tuning process

## Series: bf_train 
## ARIMA(0,1,0) with drift 
## 
## Coefficients:
##        drift
##       0.0062
## s.e.  0.0029
## 
## sigma^2 = 0.003695:  log likelihood = 616.63
## AIC=-1229.26   AICc=-1229.24   BIC=-1221.06
## 
## Training set error measures:
##                        ME       RMSE        MAE         MPE     MAPE      MASE
## Training set 2.872101e-06 0.06064977 0.03945284 -0.09724983 1.775708 0.2673914
##                     ACF1
## Training set -0.02128669

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(0,1,0) with drift
## Q* = 13.479, df = 23, p-value = 0.9408
## 
## Model df: 1.   Total lags used: 24
## Series: bf_train 
## ARIMA(0,1,0)(0,0,2)[12] with drift 
## 
## Coefficients:
##         sma1    sma2   drift
##       0.0464  0.1137  0.0062
## s.e.  0.0558  0.0576  0.0033
## 
## sigma^2 = 0.003672:  log likelihood = 618.84
## AIC=-1229.68   AICc=-1229.59   BIC=-1213.28
## 
## Training set error measures:
##                        ME       RMSE       MAE         MPE     MAPE      MASE
## Training set 4.510686e-05 0.06032785 0.0392854 -0.08257762 1.766861 0.2662566
##                    ACF1
## Training set -0.0103877

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(0,1,0)(0,0,2)[12] with drift
## Q* = 10.496, df = 21, p-value = 0.9717
## 
## Model df: 3.   Total lags used: 24

With its default parameters, the automated tuning process fitted an ARIMA(0,1,0) model that includes a drift term. With a robust search, the automated tuning process fit a ARIMA(0,1,0)(0,0,2)[12] with drift. Performance of both models are nearly identical in terms of error metrics and residual behavior. Therefore, ARIMA(0,1,0) with drift is the parsimonious model, having an MAPE score of 1.78%, and residual analysis indicating that the random component is white noise.

Linear regression with ARIMA errors

We will next train a linear regression model having the following three predictors: trend, 12-month seasonal lag, and a categorical variable for month of the year. Additionally, the errors will be modeled using the ARIMA procedure.

## Series: bf_train 
## Regression with ARIMA(1,0,0) errors 
## 
## Coefficients:
##          ar1  intercept  monthFeb  monthMar  monthApr  monthMay  monthJun
##       0.9879     0.8381    0.0029   -0.0021   -0.0011   -0.0027    0.0104
## s.e.  0.0069     0.3521    0.0095    0.0128    0.0150    0.0164    0.0171
##       monthJul  monthAug  monthSep  monthOct  monthNov  monthDec        
##        -0.0095   -0.0017   -0.0085   -0.0220   -0.0061   -0.0204  0.0064
## s.e.    0.0174    0.0171    0.0164    0.0151    0.0130    0.0097  0.0012
##             
##       0.0166
## s.e.  0.0591
## 
## sigma^2 = 0.003617:  log likelihood = 605.22
## AIC=-1178.44   AICc=-1177.14   BIC=-1113.23
## 
## Training set error measures:
##                         ME       RMSE      MAE        MPE     MAPE      MASE
## Training set -0.0009975678 0.05993264 0.039497 -0.1711719 1.777347 0.2676907
##                     ACF1
## Training set 0.009377057

## 
##  Ljung-Box test
## 
## data:  Residuals from Regression with ARIMA(1,0,0) errors
## Q* = 14.751, df = 9, p-value = 0.09801
## 
## Model df: 15.   Total lags used: 24

The linear regression model with ARIMA errors has an MAPE score of 1.78%, and the correlation analysis indicates that the random component is white noise.

Forecasting with Machine Learning Models

We will next use several machine learning regression models. Given that the purpose is to obtain a short-term forecast of 12 months, using the entire series may add noise to the model from previous cycles. Therefore, we will instead subset the series to only include the most recent cycle, beginning in January 2010, following the end of the 2008 economic crisis.

## [1] "data.frame"
## [1] 147   2
## 'data.frame':    147 obs. of  2 variables:
##  $ ds: Date, format: "2010-01-01" "2010-02-01" ...
##  $ y : num  2.28 2.28 2.24 2.36 2.31 ...
##        ds                   y        
##  Min.   :2010-01-01   Min.   :2.240  
##  1st Qu.:2013-01-16   1st Qu.:3.289  
##  Median :2016-02-01   Median :3.725  
##  Mean   :2016-01-31   Mean   :3.605  
##  3rd Qu.:2019-02-15   3rd Qu.:4.011  
##  Max.   :2022-03-01   Max.   :4.757
##           ds     y
## 1 2010-01-01 2.279
## 2 2010-02-01 2.277
## 3 2010-03-01 2.240
## 4 2010-04-01 2.364
## 5 2010-05-01 2.309
##             ds     y
## 143 2021-11-01 4.716
## 144 2021-12-01 4.604
## 145 2022-01-01 4.554
## 146 2022-02-01 4.630
## 147 2022-03-01 4.757
## [1] "ds" "y"
##         date     y
## 1 2010-01-01 2.279
## 2 2010-02-01 2.277
## 3 2010-03-01 2.240
## 4 2010-04-01 2.364
## 5 2010-05-01 2.309

Feature engineering

We will next create new features to be used as informative input for the model.

## 'data.frame':    135 obs. of  6 variables:
##  $ date     : Date, format: "2011-01-01" "2011-02-01" ...
##  $ y        : num  2.53 2.66 2.71 2.72 2.69 ...
##  $ month    : Factor w/ 12 levels "Jan","Feb","Mar",..: 1 2 3 4 5 6 7 8 9 10 ...
##  $ lag12    : num  2.28 2.28 2.24 2.36 2.31 ...
##  $ trend    : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ trend_sqr: num  1 4 9 16 25 36 49 64 81 100 ...
##         date     y month lag12 trend trend_sqr
## 1 2011-01-01 2.533   Jan 2.279     1         1
## 2 2011-02-01 2.659   Feb 2.277     2         4
## 3 2011-03-01 2.715   Mar 2.240     3         9
## 4 2011-04-01 2.722   Apr 2.364     4        16
## 5 2011-05-01 2.694   May 2.309     5        25
## 6 2011-06-01 2.774   Jun 2.400     6        36

Training, testing, and model evaluation

To allow for model comparison, follow the same procedure used for the previous models to create the training and testing partitions. It will also be necessary to create inputs for the forecast itself.

## 'data.frame':    12 obs. of  5 variables:
##  $ date     : Date, format: "2022-04-01" "2022-05-01" ...
##  $ trend    : num  136 137 138 139 140 141 142 143 144 145 ...
##  $ trend_sqr: num  18496 18769 19044 19321 19600 ...
##  $ month    : Factor w/ 12 levels "Jan","Feb","Mar",..: 4 5 6 7 8 9 10 11 12 1 ...
##  $ lag12    : num  4.1 4.1 4.36 4.39 4.47 ...
##         date trend trend_sqr month lag12
## 1 2022-04-01   136     18496   Apr 4.096
## 2 2022-05-01   137     18769   May 4.101
## 3 2022-06-01   138     19044   Jun 4.357
## 4 2022-07-01   139     19321   Jul 4.388
## 5 2022-08-01   140     19600   Aug 4.468
## 6 2022-09-01   141     19881   Sep 4.504

Model benchmark

We will use a linear regression model as a benchmark for the machine learning models.

## 
## Call:
## lm(formula = y ~ month + lag12 + trend + trend_sqr, data = train_df_st)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.41738 -0.13046 -0.05507  0.10412  0.75341 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.702e+00  2.580e-01   6.597 1.60e-09 ***
## monthFeb     1.448e-03  1.079e-01   0.013    0.989    
## monthMar     1.484e-02  1.079e-01   0.138    0.891    
## monthApr     1.171e-02  1.107e-01   0.106    0.916    
## monthMay     4.501e-02  1.109e-01   0.406    0.686    
## monthJun     8.818e-02  1.107e-01   0.797    0.427    
## monthJul     4.707e-02  1.107e-01   0.425    0.671    
## monthAug     3.269e-02  1.107e-01   0.295    0.768    
## monthSep     3.185e-02  1.107e-01   0.288    0.774    
## monthOct     1.862e-02  1.108e-01   0.168    0.867    
## monthNov     2.655e-02  1.107e-01   0.240    0.811    
## monthDec     7.783e-03  1.109e-01   0.070    0.944    
## lag12        4.609e-01  1.123e-01   4.102 7.96e-05 ***
## trend        8.244e-03  5.450e-03   1.513    0.133    
## trend_sqr   -3.874e-05  3.462e-05  -1.119    0.266    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.253 on 108 degrees of freedom
## Multiple R-squared:   0.71,  Adjusted R-squared:  0.6724 
## F-statistic: 18.89 on 14 and 108 DF,  p-value: < 2.2e-16
## [1] 0.1004068

## 
##  Breusch-Godfrey test for serial correlation of order up to 18
## 
## data:  Residuals
## LM test = 108.75, df = 18, p-value = 5.372e-15

The residual plot shows a nonrandom pattern. The correlation plot shows that the series is dependent on its lags. The density plot shows a right-skewed distribution. Therefore, it is not a valid forecasting model. However, we will use its MAPE score of 10.04% as a benchmark for the performance of the machine learning models.

Starting an h2o cluster

##  Connection successful!
## 
## R is connected to the H2O cluster: 
##     H2O cluster uptime:         1 hours 4 minutes 
##     H2O cluster timezone:       America/Los_Angeles 
##     H2O data parsing timezone:  UTC 
##     H2O cluster version:        3.36.0.4 
##     H2O cluster version age:    1 month and 15 days  
##     H2O cluster name:           H2O_started_from_R_keoka_jco123 
##     H2O cluster total nodes:    1 
##     H2O cluster total memory:   0.77 GB 
##     H2O cluster total cores:    4 
##     H2O cluster allowed cores:  4 
##     H2O cluster healthy:        TRUE 
##     H2O Connection ip:          localhost 
##     H2O Connection port:        54321 
##     H2O Connection proxy:       NA 
##     H2O Internal Security:      FALSE 
##     R Version:                  R version 4.1.2 (2021-11-01)
## 
  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |======================================================================| 100%
## 
  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |======================================================================| 100%
## 
  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |======================================================================| 100%

Now that the data has been loaded into the working cluster, we can begin the training process.

Forecasting with the Random Forest model

Build a forecasting model with the Random Forest (RF) algorithm.

## 
  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |==============                                                        |  20%
  |                                                                            
  |======================================================================| 100%
## Model Summary: 
##   number_of_trees number_of_internal_trees model_size_in_bytes min_depth
## 1              25                       25               23588        10
##   max_depth mean_depth min_leaves max_leaves mean_leaves
## 1        15   11.80000         59         84    70.60000

## 
  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |======================================================================| 100%
## [1] 0.1148354

The Random Forest model with its default settings has an MAPE rate of 11.48%, which is higher than the 10.04% error rate of the benchmark model.

Forecasting with the GBM model

## 
  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |=                                                                     |   1%
  |                                                                            
  |==========================                                            |  37%
  |                                                                            
  |==============================================                        |  66%
  |                                                                            
  |===================================================                   |  72%
  |                                                                            
  |=====================================================                 |  76%
  |                                                                            
  |=======================================================               |  79%
  |                                                                            
  |=========================================================             |  82%
  |                                                                            
  |============================================================          |  86%
  |                                                                            
  |================================================================      |  92%
  |                                                                            
  |===================================================================   |  95%
  |                                                                            
  |===================================================================== |  98%
  |                                                                            
  |======================================================================| 100%
## Model Summary: 
##   number_of_trees number_of_internal_trees model_size_in_bytes min_depth
## 1              25                       25               23588        10
##   max_depth mean_depth min_leaves max_leaves mean_leaves
## 1        15   11.80000         59         84    70.60000

## 
  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |======================================================================| 100%
## [1] 0.1017844

The Gradient Boosting Machine model with its default settings has an MAPE rate of 10.18%, which is higher than the 10.04% error rate of the benchmark model.

Prediction for March 1, 2023

## Linear regression models
data.frame(forecast(bf_lr1,h=24))[24,] # 4.819455
##          Point.Forecast    Lo.80    Hi.80    Lo.95    Hi.95
## Mar 2023       4.819455 4.491481 5.147429 4.317224 5.321687
## Exponential smoothing models
data.frame(forecast(bf_hw1,h=24))[24,] # 4.228599
##          Point.Forecast    Lo.80    Hi.80    Lo.95    Hi.95
## Mar 2023       4.228599 3.819427 4.637771 3.602825 4.854373
data.frame(forecast(bf_hw2,h=24))[24,] # 3.827890
##          Point.Forecast    Lo.80    Hi.80    Lo.95    Hi.95
## Mar 2023        3.82789 3.394173 4.261607 3.164577 4.491203
## ARIMA models
data.frame(forecast(arima_m1,h=12))[12,] # 1.557115
##          Point.Forecast   Lo.80    Hi.80    Lo.95    Hi.95
## Mar 2023       1.557115 1.46019 1.654039 1.408882 1.705348
data.frame(forecast(arima_m2,h=12))[12,] # 1.556631
##          Point.Forecast    Lo.80    Hi.80    Lo.95    Hi.95
## Mar 2023       1.556631 1.460782 1.652479 1.410043 1.703218
data.frame(forecast(arima_m3,h=12))[12,] # 1.556562
##          Point.Forecast    Lo.80   Hi.80    Lo.95   Hi.95
## Mar 2023       1.556562 1.460844 1.65228 1.410174 1.70295
data.frame(forecast(arima_a1,h=24))[24,] # 4.190090
##          Point.Forecast    Lo.80    Hi.80    Lo.95    Hi.95
## Mar 2023        4.19009 3.808458 4.571722 3.606434 4.773745
data.frame(forecast(arima_a2,h=24))[24,] # 4.211082
##          Point.Forecast    Lo.80    Hi.80    Lo.95    Hi.95
## Mar 2023       4.211082 3.821692 4.600473 3.615561 4.806603